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Abstract 

The density density correlation function is computed for the BogoHubov pseudoparticles created 
in a Bose-Einstein condensate undergoing a black hole flow. On the basis of the gravitational 
analogy, the method used relies only on quantum field theory in curved spacetime techniques. 
A comparison with the results obtained by ab initio full condensed matter calculations is given, 
confirming the validity of the approximation used provided the profile of the flow varies smoothly 
on scales compared to the condensate healing length. 
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I. INTRODUCTION 



It is by now well known that the behavior of elementary excitations in various condensed 
matter systems can be described, under some approximation, by an effective theory of mass- 
less scalars (phonons) propagating on a fictitious curved spacetime according to a covariant 
Klein-Gordon equation [HE]- This curved spacetime, which has nothing to do with the real 
physical spacetime on which the underlying condensed matter theory is defined, is deter- 
mined by the background features of the condensed matter system. In particular, curvature 
is associated with inhomogeneities of the system. The approximation, which is at the core of 
this condensed-matter analogy, is typically a hydrodynamical approximation whose validity 
requires one to consider the system on scales much larger than the typical length associated 
with the microscopic (atomic) structure of the system. 

It is of particular interest to consider configurations for which the associated curved 
spacetime metric describes what in gravity would be called a black hole, BH. This happens, 
for example, in a stationary flow when the speed of sound varies with position in such a way 
that there is both a subsonic and a supersonic region with the direction of flow from the 
subsonic towards the supersonic region. The surface separating the subsonic from the super- 
sonic regime acts as a sort of "acoustic horizon" since sound waves in the supersonic region 
(the "acoustic" BH) are trapped inside this horizon being unable to propagate upstream 
towards the subsonic part, mimicking the trapping of light (and everything else) inside the 
horizon of a gravitational BH. 

There is therefore the concrete possibility of investigating in these condensed matter 
configurations the analogue of one of the more exotic processes foreseen in theoretical physics, 
namely BH quantum evaporation as predicted by Hawking in 1974 [3]. This process, despite 
what one would naively think, is not peculiar to gravity. It is completely kinematical as it 
relies only on the presence of a horizon in the underlying metric on which the quantum fields 
propagate. As such it is expected to occur in condensed matter analogue systems when the 
associated curved spacetime has a BH form [H [2] . 

Hawking evaporation is a general pair creation process in which quantum vacuum fluc- 
tuations are converted into real on shell quanta. One quantum (the positive energy one) is 
emitted outside the horizon and propagates to infinity where it is characterized by a thermal 
spectrum at the Hawking temperature, the other member of the pair (the negative Killing 
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energy quantum), called its "partner" [1], propagates inside the horizon and remains trapped 
there. It is this compensation of positive and negative energy which allows the process to 
occur in stationary, even static configurations. 

For a gravitational BH the Hawking temperature associated with this emission process 
is extremely low Th ~ 10""^ {Msun/M) K, where M is the BH mass. This makes Hawking 
evaporation astrophysically irrelevant and impossible to detect for stellar mass or larger 
black holes since the signal is completely overwhelmed by other sources including the 2.7 K 
cosmic microwave background radiation. In principle Hawking radiation from primordial 
BHs E] with small enough masses could be detected [H E], but no BHs of this type have 
been discovered. 

Thus analogue systems appear today as the only realistic hope of finding experimental 
verification of Hawking's prediction. But even in this context the experimental situation is 
difficult due to competing effects such as thermal fluctuations and quantum noise. These ef- 
fects can mask the signal even in the most favorable situations where the associated Hawking 
temperature is as high as 1/10 of the background temperature. 

A way to bypass this critical problem was proposed in [9j where it was shown that the 
Hawking process encodes characteristic correlations between the quanta and their partners 
that can be measured in analogue systems, allowing for clear identiflcation of the Hawking 
process. This kind of in-out correlation is not of physical interest for a gravitational BH 
since external observers have no access to the interior BH region (where the partners live). 
For acoustic BHs the region inside the horizon is simply the supersonic part of a flow and 
hence is physically accessible for measurements. The relationships between the equal-time 
correlations studied in analogue BHs [9l [10] and the space-time correlation pattern found in 
relativistic theories are discussed in [TT] . 

One of the most highly studied systems which can be used to create an analogue black 
hole is a Bose- Einstein condensate [12], BEG, which is an ensemble of sufficiently cold atoms, 
the vast majority of which are in the ground state of the system. The small fraction of the 
atoms which are in excited states can be described as (quantum) excitations above a classical 
(c number) condensate field. 

In Ref . [9] , using the condensed matter - gravity analogy, the correlations associated with 
Hawking radiation in a BEG elongated in one spatial direction undergoing a BH-like fiow 
were computed using the standard tools of quantum field theory, QFT, in curved space. 
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Two approximations were made to simplify the calculations. The first was to reduce the 
mode equation from a 4-D form to a 2-D form. In the dimensional reduction process there 
is a potential that appears in the 2-D mode equation which causes backscattering of the 
modes to occur. The second approximation was to neglect this potential. Without it the 
modes propagate freely as 2-D plane waves. 

It was found that if a sonic horizon is present then there is a negative correlation peak 
in the density density correlation function when one point is inside and one point is outside 
the horizon. No such peak occurs if there is no horizon. In [TU] ab initio full condensed 
matter calculations of the density-density correlation function were made. The calculations 
were based on Monte Carlo simulations within the full microscopic quantum description of 
the Bose-Einstein condensate. Comparisons were made with the analytic calculation in |9] 
and it was found that the two calculations are in approximate agreement for both the size of 
the peak and the width of the peak so long as the flow varies smoothly on scales compared 
to the condensate healing length. 

The analytic calculation in P] was carried out only for the case in which one point is 
inside and one point is outside the horizon. However, it can be generalized in a straight- 
forward way to consider all possible pairings of the two points. When this is done the only 
other structure found is the usual peak that occurs when the two points come together. 
The calculation of the density density correlation function in |10] showed a richer structure. 
In particular a negative correlation peak was found in the case that both points are inside 
the horizon. This is a second effect which depends on the existence of a horizon. A third, 
positive correlation peak was also found [ini [131 E] when one point is inside and one point 
is outside the horizon. It was found that this peak occurs even if the horizon is not present. 

Because the calculation in j9] neglected the potential in the 2-D mode equation it is 
natural to ask whether this potential has any important effects. In this paper we answer 
that question by numerically computing the density density correlation function when the 
potential is included in the mode equation. We find that the resulting backscattering of the 
modes leads to a more complex structure for the correlation pattern which fits very well 
with the one obtained in |10j . 

In Sec. II the relevant properties and equations for a BEC are reviewed. In Sec. Ill the 
specific model we consider is discussed. In Sec. IV the Unruh state which we use for our 
calculations is discussed and formulas for the relevant Bogolubov coefficients are derived. In 
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Sec. V the derivation of the specific expressions for the density density correlation function 
that are used in our numerical calculations are given. Sec. VI contains a derivation of the 
specific mode equations which we solve. In Sec. VII the numerical computation of the density 
density correlation function is discussed. The results of the numerical calculations are given 
in Sec. VIII and our conclusions are given in Sec. IX. 



II. BOSE-EINSTEIN CONDENSATES AND THE ASSOCIATED ANALOGUE 
METRIC 

Here we shall briefly review the basic equations describing a BEC and how these lead, 
under the hydrodynamical approximation, to the context of an acoustic metric effectively 
governing the propagation of the fluctuations. 

In the Bogoliubov theory the basic bosonic-field operator ^ describing the atoms is split- 
ted in a classical mean field \l/o describing the macroscopic occupied low energy state of the 
system (the condensate) and a part describing the quantum fluctuation above this classical 
state 

^ = ^o(l + 0). (2.1) 
The evolution of the condensate is governed by the Gross-Pitaevskii equation 

ihdr-^^ = f-T^^V' + Ve^t + ^7|^orl *o , (2.2) 
\ zm / 

where m is the mass of the atoms, g the coupling and Vext the external potential confining the 
atoms. The field operator describing the noncondensed part satisfies the Bogoliubov-de 
Gennes, BdG, equation 

Zm m Wo 
where n = |\I^oP is the condensate density. 



The BdG equation (2.3) can be manipulated to obtain a curved space wave equation as 
follows. One rewrites the condensate wave function \l/o and the basic field operator \l/ in the 
so called density-phase representation, namely 

if = J^^Tfh e'^^+^'^ - ^o(l + — + tOi) . (2.4) 

2n 



In this representation the BdG equation becomes a pair of equations of motion for the 
density fluctuation fii and the phase fluctuation §1 



hdT0i 



mc 



n 



-fii + 



^ e vhv(— )] = 

An ^ ^ n'^ 



-V(woni + — V^i) 
m 



(2.5) 
(2.6) 



where vq = hVO/m is the condensate flow velocity and c = \frigjrn the local speed of sound 
which will play a critical role in our construction. £^ = h/mc is the so called healing length, 
it is the fundamental scale in describing the microscopic quantum structure of the BEG. 
If we limit ourselves to scales much bigger than ^, i.e. within the so called hydrodynamical 



approximation, the last term in Eq. (|2.5|) can be neglected namely 



hn 



rii 



(2.7) 



Substituting this equation into Eq. (2.6) then decouples the system of equations (2.5, 2.6) 
yielding the equation 



-* Tl -* -* Ti -* 

- {Ot + Vvo) — -{Ot + voV)ei + V-VOi = 



m 



which can formally be rewritten as 



□^1 = 



where 



□ 



--d., 
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and g'^'^ is the inverse of the 4x4 matrix 



n 

mc 



5i 



(2.8) 

(2.9) 
(2.10) 

(2.11) 



Here g is the determinant of g^j,y and Latin indices range from 1 to 3 and are used for spatial 
components. 



Eq. (2.8) therefore has the form of the covariant Klein-Gordon equation for a massless 



scalar field propagating in a fictitious curved spacetime whose line element is 



n 



mc 



\—(?dT'^ + [dx — vodT){dx — VQdT)\ 



(2.12) 



This is the core of the analogy: within the hydrodynamical approximation, the Bogoliubov 
theory of phase fluctuations in a BEG is equivalent to QFT in curved spacetime for a mass- 



less scalar field. Looking at the metric (2.11) we see that for condensate configurations for 



which the flow turns supersonic, i.e. admitting a region where l^ol > c, the associated ana- 
logue metric describes what gravitational physicists would call a BH in Painvleve-GuUstrand 
coordinates. 

On the basis of the analogy one expects this kind of BEG configuration to show a Hawking- 
like process in the form of emission of correlated pairs of phonons on opposite sides of the 
horizon. The study of these correlations will be performed using the tools of QFT in curved 
spacetime for an idealized flow profile. 

The results will show the validity of this much simpler theoretical scheme to handle a 
complex system like a supersonic flowing BEG provided the background quantities vary on 
scales much larger than ^. 

III. THE MODEL 

Our model consists in an infinite elongated (along the x axis) condensate whose transverse 
size l± is constant and much smaller than ^. So the dynamics is frozen in the transverse 
direction and the problem becomes basically one dimensional. The flow is stationary and 
directed along x, from right to left, i.e. vq = —vqx with vq a positive constant. The sound 
speed is adjusted (see below) so that for a; > 0, it is Vq < c(x), while for a; < 0, fo > c(x). 
So X = plays the role of the horizon and the supersonic {x < 0) region is the BH region. 

We will call the region x > the R region and the region x < the L region. A Penrose 
diagram for the acoustic metric is given in Fig. [Tj We further assume the density to be 
constant. Our nontrivial configuration for the sound speed can be obtained by changing the 
coupling constant g (and hence the speed of sound) along x by means of a spatially varying 
magnetic field in the vicinity of a Feshbach resonance [12] . 

Quantization is performed expanding the hermitian quantum operator 9i in terms of cre- 
ation and annihilation operators satisfying usual bosonic commutation rules. Symbolically 



where the and their complex conjugates are a complete set of solutions to the mode 
equation 




(3.1) 



= , 



(3.2) 



with □ is given in Eq. (2.10). The choice of the set of solutions {/j} selects the vacuum 
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FIG. 1: Penrose diagram for an acoustic black hole metric, 
state for the field as di\0) = 0. The appropriate vacuum to describe Hawking radiation is 



discussed in the following section. The equation (3.2) once expanded has the same form as 



(2.8) with 9i replaced by /«, vq = —vqx with vq and n positive constants. The modes fi 
are assumed to be functions only of T and x. Our assumption regarding the scale of the 
transverse dimensions ^ ^) forbids excitations with transverse momenta. 

In order to solve the mode equation a sequence of coordinate transformations, familiar 
in BH physics, will be performed to simplify it. First we introduce a "Schwarzschild-like" 
time t as follows: 

xi 



t = T- dy , (3.3) 



in the R region, while in the L region 



The constants xi,X2 and a are arbitrary. 
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One can then write the mode equation as 



[□(2) + v)df'^ = 



where 9\ ^ is related to 9i by 



and 



□ (2) 

V 



c d'^ — vl d"^ dc 

1 u 1 

— f Q dt"^ c dx"^ dx 



1 + 



d_ 
dx 



Id^c 
2dx^ 



1 - 



Vn 



dc 



+ 



5fn / dc 



c^ J Ac \dx J ' Ac^ \ dx ^ 
Note that D*^^^ is the covariant d'Alambertian for the two dimensional metric 



(3.5) 
(3.6) 

(3.7) 
(3.8) 



ds' 



'-de + 



rdx 



C^ — Vn 



(3.9) 



Examination of (3.9) shows that unhke the coordinate T which is always timelike, the 



coordinate t is timelike in the R region and spacelike in the L region. It approaches the 
value — oo on the past horizon, H_, and the value +oo on the future horizon H^. 

It is useful to define in both the R and L regions the "tortoise" coordinate x*. In the R 
region the definition is 

civ) 



Jx3 c^[y)-v^ 

In the L region the definition is the same but with a different, in general, lower limit which 
we will call 0:4. In the R region x* ranges from —00 on the past and future horizons {H^) to 
+00 in the limit x — > +00. In the L region x* is again —00 on the past and future horizons 
and increases to +00 in the limit x ^ —00. It thus acts as a typical time coordinate, as can 
be seen rewriting the metric as 



ds' 



c^ — v} 



^i-dt^ + dx*^) . 



(3.11) 



The utility of this tortoise coordinate is that the wave equation for the modes takes the form 



92 



92 



dt"^ dx* 2 



cff 



with 



eff 



(3.12) 



(3.13) 
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With these definitions we can define in both the R and L regions the retarded and advanced 
null coordinates respectively 

u = t — X , 

V = t + x* . (3.14) 

Because the mode functions propagate across the future horizon /"•", it is useful to have v 
be a continuous variable across that horizon. As shown in Sec. VIII one way to accomplish 
this is to choose xi = X2, X3 = X4, and then find the appropriate value for a. 

For later use we introduce another set of retarded and advanced coordinates called Kruskal 
coordinates which are defined as 

Uk = -e"''"/K, (3.15a) 
Vk = e^^K , (3.15b) 

in R region, while 

Uk = e-^^'/K , (3.16a) 
Vk = e'^Vft , (3.16b) 

in L region. Therefore Uk = on the future horizon and Vk = on the past horizon 
H~ . The parameter k is the surface gravity of the horizon, in our acoustic setting it is 

dc 

. = -Uo. (3.17) 



In Ref. [9], the term containing Ve/f in the mode equation (3.12) was neglected with the 
result that the mode solutions can be written in terms of the simple plane waves (e~™" 
and e~™^) that propagate freely without being scattered. Here the full equation (with Vg// 
included) will be solved (numerically) to also take into account the backscattering which in 
the 4-D case is caused by the effective curved geometry (i.e. the inhomogeneities of the BEC 
medium) and in the 2-D case is caused by the effective potential. 



IV. THE UNRUH STATE 



To proceed further one has to select the quantum state for the field encoding the 
formation of the acoustic BH which is the process which triggers Hawking radiation. 

11 



However if one limits our analysis to sufficiently late times after the BH formation, the 
features of the emitted Hawking quanta are completely independent on the details of the for- 
mation process, being determined only by the properties of the resulting stationary horizon, 
more specifically its surface gravity k. So one can avoid discussing the complicated physics 
underlying the dynamics of the formation process and simply impose on the past horizon 
H~ of our stationary BH metric the appropriate, but universal, boundary conditions which 
mimic the horizon formation process. 

This is the spirit of the so called Unruh state [T3]. Technically Hawking radiation is 
superimposed on our stationary background by requiring the retarded modes originating 
from the past horizon H to be positive frequency with respect to the Kruskal coordinate 
Uk which is an affine parameter along H~ . 

For this state there is a flux of radiation which comes through past horizon, H~ while at 
past null infinity, i.e. there is no radiation. To obtain this behavior the scalar field 9{ 
is expanded in terms of two sets of modes so that 



K(, , ^\ t J I, ,\„,K*t 



du * {u,x) + * *{oo,x)] . (4.1) 

Here both sets of creation and annihilation operators obey the usual commutation relations 
and both sets of modes are solutions to the mode equation 

+ ^ + Vos]uiu,x) = 0. (4.2) 



The modes uf originate at I]^. They are positive frequency with respect to the time coordi 



nate t and on have the form Uj ~ e^^'^'". Because of the potential term in Eq. (4.2) they 
are partly transmitted towards and partially reflected to and see Fig. [2] 

The modes which come through H are, on , positive frequency with respect to 
the Kruskal null coordinate Uk- There they have the form ~ e"**^^^^. They are partially 
transmitted to or and partially reflected to I^, see Fig. [s} 

Note that these behaviors on and H~ define the state because it is always possible 
to fix the mode function and its first derivative in any way that one wishes on a Cauchy 
surface. Different ways of fixing these initial conditions lead to different states for the field. 

The modes are normalized using the scalar product [16] 

01 (a;) (l);{x)^/-g^{x)n>'dE , (4.3) 
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Modes 




FIG. 2: Incoming modes uf from 

with E a Cauchy surface and a future directed unit vector which is perpendicular to E. 
The normahzation condition is 

{u{uj,x),u{uj\x)) = 5{uj — oj') . (4.4) 

Normahzation of the ui modes is done on the surface . Technically one needs a Cauchy 
surface. A full Cauchy surface for the R region would consist of along with the part of 
H which bounds this region. Then one could add to this a spacelike or null surface covering 
the rest of the spacetime. What is chosen is unimportant because the modes which originate 
on are zero on the past horizon and along the above mentioned surface. On the surface 
the normalization condition is 

dvuT{uj,x) dv uT*{i^'.x) = 5{uj-u') . (4.5) 

-oo 

The result is that on I^, 

ur\uj,x)^^=. (4.6) 
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Modes u , 




FIG. 3: Outgoing modes from H . 

Normalization of the the modes is done along H~ which is a Cauchy surface. The 
normalization condition is 



— oo 



The result on H is 



(4.8) 



Although one can define initial data for the modes on H~ , it is difficult to evolve these 
modes because in the Kruskal coordinates Uk and Vk the mode equation is not separable 
and in the t and x coordinates it is not easy to express the initial conditions for these modes. 
The mode solutions discussed above that begin on I]^ are zero on H . However, there exists 



a set of solutions to the mode equation (4.2) which are positive frequency with respect to 



t and which pass through the part of H~ which intersects the R region. On that horizon 
they have the behavior ~ e"""^". 

These solutions are zero on both and the part of H~ which borders the in region. 



14 



They are normalized on H using the condition 

{uHiuj,x),ufj{uj',x)) = -i duu^{u,x) duufj*{u}',x) = 5{uj - cu') , (4.9) 



with the result that on H 



Mg(u;,x) = -=, (4.10) 
v47ra; 



see 



Fig-i 




FIG. 4: Ougoing modes u^, defined in the exterior region [x > 0). 

For the part of H~ which serves as a boundary for the L region one has a complete set 
of modes which on H~ are positive frequency with respect to the time coordinate x*. They 
are zero on the part of H~ which borders the R region. On the other part of they have 



the form ~ e*'^" and are normalized using the condition (4.9) with the result that on 
the part of H~ which borders the in region, 



<(a;,x) = -=, (4.11) 



see Fig. [5] 
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Modes 




FIG. 5: Outgoing modes u^, defined in the interior region {x < 0). 

All that remains is to express the modes in terms of the modes and u^. This 
is done through Bogolubov transformations whose coefficients can be computed using the 
scalar product evaluated on the past horizon H~ . We can write 

POO 

Jo 

poo 

+ ■ (4-12) 

Jo 

The coefficients can be obtained by using the scalar product evaluated on H~. RecaUing 
that u^- = on the part of H" which borders the L region and u^- = on the part of H~ 
which borders the R region we have 

«<?K0. = {uh{^^k,x),u^{u},x)) , (4.13a) 

= -{uU^K,x),u^H*ico,x)) , (4.13b) 

(^^K^ ^ ('^h{^k,x),Uh{uj,x)) , (4.13c) 

= -iu^iujK,x),u^H*i^,x)) ■ (4.13d) 
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Using the expressions (4.8) and (4.10) and evaluating the scalar product on H one finds 

. (4.14a) 



a. 



R 



(HIk exp 



-iukUk - i-\og[-UK) 

K 



lUJ 



+ IUJk 



K) 



Making the change of variables Z = —Uk allows one to compute the integral in terms of 
Gamma functions with the result that 



a 



R 



1 



2'n-K Y lok 



(4.14b) 



In a similar way one finds 



R 



.UJ 



(HIk exp —iojKUK + i—^og(—UK) 

L K 



a 



LdK ^ 



(IUk exp —iukUk + log(f//^ 



-lUJ 



+ iujK 



Pi 



2ttk 

4:71 ^ukuj Jo 
-^./^r(.a;/«:)(.u;^)-^-/^ 
i 

Zttk y Uk 



(4.14c) 



kU. 



K 



(4.14d) 



K 



-lUJ 



kUk 



(4.14e) 



V. DENSITY CORRELATIONS 



Correlations are a manifestation of the entanglement existing between created particles 
and their partners [9l [10] , see also [T71 - [T9j . As discussed in Sec. IV, the Unruh vacuum we 



have constructed, i.e. the state annihilated by the operators ax and h in Eq. (4.1) describes 



Hawking radiation at late times. The created quanta appear at late times as the three types 
of modes which are depicted in Figs. |6l [7l and [8] 



The first mode (Fig. 6) at late times for x — )• oo has form and describes phonons 
propagating towards x — )■ +oo with velocity V = (c+ — vq) > where c+ = c(x = +oo). 
The second mode (Fig. 7) at late times for a; — t- — oo has the form and describes 
phonons propagating inside the BH towards x — )■ — oo with velocity V = (c_ — vq) < 
where c_ = c{x = — oo). These are the so called partners (negative Killing frequency). The 
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FIG. 6: Phonons (Hawking quanta) created in the exterior (subsonic) region. 




FIG. 7: Negative frequency phonons (partners) created in the interior (supersonic) region. 

last one (Fig. [s]) also describes modes propagating inside the BH, this time with velocity 
V = — (c_ + t>o) < 0; their asymptotic {x — — oo) form is ^^7=- (positive frequency). 

Two-point correlations related to Hawking radiation are between these modes. There 
are therefore three kinds of relevant correlations. The first is L — R and correlates quanta 
between the modes in Fig. |6] and Fig. [7], i.e. modes on opposite sides of the horizon, see 
Fig. [9} Associated with this one expects to observe a peak for the correlations in the {x, x') 
plane along the line 

= (5.1) 

C_ -Vq C+ - Vq 

with x' < and x > 0. There is another correlation of the L — R form which correlates 



quanta between the modes in Fig. [6] and Fig. [8| see Fig. 10 The peak is along 
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FIG. 8: Positive frequency phonons created in the interior (subsonic) region. 




FIG. 9: Left-Right correlations between phonons in Figs. [6| [7[ 



^ ^ (5.2) 



Finally there is a L — L correlation (both modes inside the horizon) of the modes in Fig. [T] 



and Fig. [8| see Fig. [TTJ The peak this time is along 

(5.3) 



x' X 



-(c_+fo) C+-VO 

where both x, x' < 0. 

The equal time density density correlation function, which is the basic object of our 
investigation, is formally defined as 

G2(T, x; T', x') = lim (ni(T, x)ni{T', x')) . (5.4) 
T->r' 
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FIG. 10: Left-Right correlations between phonons in Figs. [6j [8} 




FIG. 11: Left-Left correlations between phonons in Figs. [7j [8j 



Writing fii in terms of the phase operator Eq. (2.7) one can relate G2 to the symmetric 

^ (2) 

two-point function of the field 0^^ 

hn 



G2iT,x;T',x' 



2m£j_c^(x)c^(a;') t'-s-t 



- lim D^c{x)c{x'){{d'i'\t,x),e'i'\t',x')}) , (5.5a) 



(5.5b) 



where {, } is the anticommutator. Here one should use the definition (3.3) for t in the R 
region and the comparable definition in the L region. 
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The two point function is computed using Eqs. (4.1), with the result that 

{{e'i'\t,x),9?\t',x')}) = i + j, 



(5.6a) 



J 



poo 

/ duK [uH{0JK,t,x)uH*{u}K,t',x') +UH*{uJK,t,x)uH{uJK,t',x')] , (5.6b) 

Jo 

du [uf{u,t,x)uf*{uj,t',x') +uf*{u,t,x)uf{u,t',x')] . (5.6c) 



Substituting Eq. (4.12) into Eq. ( |5.6b ) and using Eqs. ( |4.14[ ) one finds 16 integrals of the 
form 

i / , I.I. i\UJUJ 



f^C<0 POO POO 

, duK / dco duj 
An^K^ .In .In .In OJR 



JO JO 

±lUJK)^''''''{±iuJK)^'''''''Ui{uj, t, X)U2{0J', t', X') 



(5.7) 



with all combinations of ± and Ui and U2 being various combinations of u'j^, u^, and their 
complex conjugates. The integral over uk can be computed. With the change of variable 
z = (Inwii-)/^ the result is proportional to either 6{ijJ — u') or 6{ijJ + ijj'). Given the limits of 
integration terms containing 6{uj + u') integrate to zero. Using the relation 

-iU)\ 7CK 



\ K, J \ K 



one finds that 
/ = 



cjsinh (^) 

[u'^{u, t, x) U^{u, t', x') + M^*(w, t, x) U^{(jJ, t', x') 



(5.8) 



OO 

duj 7 ^ 

sinh (^) 

+U^{uj, t, x) U^{bJ, t', x') + U^*{u, t, x) U^{uj, t', x') 

+ cosh (^~^ [^ni^^ '^^ ^) "^H ('^) + "^iii^j ^) ^) '^ni^^ 
+ii|(w, t, x) u^*{u, t', x) + uf{uj, t, x) u^{u, t', x')] } 



(5.9) 



VI. MODE FUNCTIONS 



It is next necessary to obtain explicit expressions for the mode functions. Their normal- 
ization has been described above. For each of them the mode equation can be solved using 
separation of variables. The result is 

1 



UH{uj,t,x) 
UHiuJ,t,x) 



uf{uj, t, x) 



y/Antu 
1 



1 



\/ 4:700 
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(6.1a) 
(6.1b) 
(6.1c) 



In all cases the part of the mode function that depends on the coordinate x obeys the 
equation 



+ io'x + V,sX = 0. (6.2) 



This equation must be solved numerically in general. If c — )■ constant as x — > ±00, then 
Ves — )• in these hmits. Further, because c — f at the horizon (x = 0), Vcs — )■ at the 
horizon as well. This aids in determining the boundary conditions on x- 

The easiest case to consider is Xh- region, x* is the time coordinate. Thus on 

both the past and future horizons (which are both at x = 0) one has x = e"*"^^*. Since Xh 
enters the L region from the past horizon H~ the modes are initially right moving. This is 
the reason for the factor of e*'^* in Eq. (6.1a). For numerical purposes it is useful to break 



Xh i^to real and imaginary parts. Thus we define Xc be the solution to Eq. (6.2) which 
in the limit x — i- 0^ has the behavior 



Similarly in this limit we define 



Then 



Xc cos{ux*) . (6.3a) 
Xs -> sm{ux*) . (6.3b) 
Xh = X'c + . (6.3c) 



The numerical computation of Xc Xs described in more detail in the next section. 
In the R region the situation is more complicated. It is useful to define two solutions to 



Eq. (6.2) by their behavior in the large x limit: 

xf — )■ cos(a;x*) , (6.4a) 
xf -> sin(cux*) . (6.4b) 

Near x = these solutions will have the form 

Xc Acos{uJX*) + B sm{ujx*) , (6.5a) 
xf — Ccos(a;x*) + sin(a;x*) . (6.5b) 

Using these, one can define solutions which correspond to right moving and left moving 
waves in the large x limit as 

xT = Xc+^xf, (6.6a) 

xT = Xf - ^Xf ■ (6.6b) 
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Near x = these solutions will have both right and left moving parts so that 

Ere'^""* + F^e"^"^^* , (6.7a) 
xT ^^e'""* + i^^e"'""* . (6.7b) 



One easily finds that 



Er = ^[A + D-i{B-C)], (6.8a) 

Fr = ^[A-D + i{B + C)] , (6.8b) 

Ee = ^[A-D-t{B + C)], (6.8c) 

Fe = ^[A + D + i{B-C)]. (6.8d) 



For the modes which enter the R region from the past horizon, the normalization occurs 
at that horizon while the boundary condition on Xh i^sX the mode function should be a 
right moving wave in the limit x — i- oo. Thus 



Xh 



Nx 



oo 
r 1 



(6.9) 



The normalization constant is determined through the normalization condition (4.10) on 



H . Near x = it is the right moving part of Xh which corresponds to the part of the mode 



function coming from H . Using Eqs. (6.1b) and (6.7a) one finds that 



N 



Er 



(6.10) 



To find the behavior of for x < we note that it is the left moving part near x = O"*" 
which goes through the future horizon. Then by continuity of the total mode function 
across the future horizon one finds that for x < 



Xh 



Er 



iXc-^Xs) ■ 



(6.11) 



Finally we discuss the modes uf which originate on / in the R region. The normalization 



for these modes occurs on / and is given in Eq. (4.6). Because of backscattering part of 
the mode function reaches I^. Therefore, 



xf = xT + Kx 



oo 
r ! 



(6.12) 
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with K a constant. The boundary condition is that on H , uf = 0. This is accomphshed 
by having Xh ^'g a left moving wave in the hmit x 0"^, that is 

xf = Ge-'^^* , (6.13) 



with G another constant. Using Eqs. (6.7) in Eq. (6.12) and setting the result equal to (6.13) 
gives 

K = (6.14a) 
G = Fl-^. (6.14b) 
For X < we again use continuity of the mode function at to obtain 

This ends the derivation of the equations necessary for the computer program. We used 
Mathematica to combine them to give expressions for the density-density correlation function 
in terms of the modes Xc-> X^i xfi cind their derivatives and also in terms of the constants 
A, B, C, D. All of these quantities are numerically computed. 



VII. NUMERICAL COMPUTATIONS 



In this section we discuss the numerical computation on the density density correlation 
function. This involves numerical computations of the mode functions in both the L and R 
regions. The results are substituted into the expressions for the density density correlation 
function which involve integrals over the frequency u. Because of the different behaviors 
of the mode functions in the in and out regions there are different specific expressions for 
the density density correlation function in terms of the mode functions and their derivatives 
depending on where the two points are located. These have been obtained using the algebraic 
manipulation program Mathematica but they are too long to be shown here. 

There is a problem that arises in computing the mode integrals in all but one case. That is 
that the process of computing the integrals over uj and computing the derivatives of the two 
point function necessary to obtain the density density correlation function do not commute. 
The point is that for the two point function the integrals over u are well defined and converge 
in the limit a; — )■ oo so long as the points are separated. However, if one takes the relevant 
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derivatives first and thus obtains an integral over derivatives of tlie mode functions, then the 



resulting integrals over u are not well defined. This is discussed in more detail in Sec. VII B 



where it is shown that a subtraction procedure can be used to overcome this problem. Note 
that there is also the usual infrared divergence of the two point function in two dimensions, 
but this can be taken care of with a judicious choice of infrared cutoff, which in itself is 
reasonable because of the finite dimensions of the physical system being modeled. 



A. Numerical Calculations of the Radial Mode Functions 



In both the L and R regions the numerical computation of the radial mode functions can 
be accomplished by computing the functions Xc and Xs for these regions. In the R region 
these mode functions asymptotically approach cosux* and sincjx* as a; — )■ oo. In the L 
region they have the same asymptotic behaviors in the limit x — )■ 0~. The starting values 
for these modes can be made arbitrarily accurate by starting at a large enough value of x 
for the R modes and a value of x close enough to zero for the L modes. 

The only other task that needs to be accomplished before the mode functions are substi- 
tuted into the expressions for the density density correlation function is that the parameters 
A, B, C, and D must be numerically determined. This is accomplished by matching the 



numerically computed mode functions Xc X7 to Eqs. (6.5). One way to do this is to do 
the matching at a number of values of x with < a; ^ 1 and then use a linear extrapolation 
routine to find the values that A, B, C, and D in the limit that x — > 0. 



B. Numerical Computation of the Mode Integrals 

The computation of the mode integrals in the density density correlation function would 
be a straight-forward numerical exercise if they converged. However there are two problems. 
One is that because we are working with an effective relativistic two dimensional quantum 
field theory there is an infrared divergence in the two point function and therefore there are 
infrared divergences in the density density correlation function. The other is related to the 
fact mentioned previously that strictly speaking one must compute the mode integral in the 
two point function and then take the required derivatives of it to obtain the density density 
correlation function. 
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Infrared divergences 



The simplest way to numerically remove the infrared divergence is to impose a lower limit 



cutoff on the mode integrals. From Eqs. (5.6a), (5.6c), (5.9), (6.1), and (6.2) one can see 
that the two point function has infrared divergences which go like both l/Xi and \og{\i) 



for a small enough infrared cutoff A^. As shown in Eq. (5.5) the density density correlation 
function is composed of terms containing two derivatives of the two point function. Each 
time derivative brings down a factor of u while each space derivative results in one term in 
which a factor of u is brought down, one term in which one or more derivatives of the sound 
speed c(x) are found, and one term in which a derivative of the spatial part of one of the 
mode functions occurs. If the sound speed at the point where the spatial derivative is taken 
is constant then in all surviving terms a factor of u is brought down. If the sound speed is 
not constant then there are one or more terms in which either no factor of u or only one 
factor of uj are brought down by the spatial derivatives. Thus if both points of the density 
density correlation function are in regions where the sound speed is constant (and therefore 
the effective potential is zero) then there are no infrared divergences in the density density 
correlation function. Otherwise infrared divergences of the types discussed above occur. 

Even if there are no infrared divergences, the finite physical size of the system means that 
there will be an infrared cutoff. As discussed above for low enough values of the cutoff there 
will be terms containing the cutoff which go like 1/A^ and log(A^). However if a higher value 
is given to the cutoff then oscillatory behavior in the density density correlation function is 
generated by the cutoff. To see why this occurs one can look at the behaviors of the mode 
functions when the effective potential is zero. Then it is not hard to show that there will 
be terms which go like ^(A^Am), ci{\i^v), and so forth. If one point is outside the horizon 
near one end of the BEG and the other is inside the horizon near the other end of the BEG 
then roughly speaking Au ~ Af ~ Ax, with Ax the size of the BEG. Since the sound speed 
is of order unity in the models we consider it is clear that in this case the arguments of the 
cosine integral functions are of order unity. This leads to oscillatory behavior in terms of 
the dependence on the value of the cutoff A^. However, the approximations we are using 
assume a BEG of infinite length and so are not expected to be valid near the edges of the 
BEG where other effects that are being neglected may be important. Therefore we shall 
restrict consideration to much smaller values of Ax such that A^|Ax| ^ 1. The effects of the 
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infrared cutoff on the density density correlation function when this restriction is satisfied 
are discussed further in Sec. IVIIII 



Large frequency behavior 

When computing the density density correlation function numerically one must first find 
numerical values for the mode functions and their derivatives and then substitute these 
into the appropriate expression for the correlation function and numerically compute the 
integrals. However, if this is done one finds that some of the resulting integrals do not 
converge. Instead for large values of the frequency u their integrands oscillate in u but with 
either growing or approximately constant amplitudes. An illustration of what is happening 
can be obtained from the cosine integral function. 

ci(Aa;) = - / du . (7.1) 

If one computes the derivative with respect to x of both sides and interchanges the order of 
integration and differentiation on the right hand side then 

cos(Ax) 



oo 



du sin{xuj) . (7.2) 



A 



X 

The right hand side of this equation is clearly not well defined since the amplitude of the 
integrand is constant in the limit a; — t- oo. 

One way around this problem is to evaluate the integrands for the integrals in the density 
density correlation function in the large u limit. Then the terms that are poorly behaved 
at large u can be identified and subtracted off. Because it is the large u limit, this can be 
done analytically and then the specific terms that are subtracted off can be added back and 
treated analytically. A simple example of how this works would be 



oo 



A 



, cos(ux) f°° , / cosiux) cos(uJx)\ , , ^, 

du , ' = du { , ^ ' ^ — - - ci(Ax) . (7.3) 



Then if one wanted the first derivative of this with respect to x one would have 

d /"°° , cos(ci;a;) /"°° , / sin(a;x) , A cos(Aa;) 
— / du , ' = - du { u , ' - sm(a;x) H ^ — - . (7.4) 

In this case the integral on the right hand side has an integrand that oscillates but with an 
amplitude that vanishes in the limit a; — )■ oo. It can therefore be computed numerically. 
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To evaluate the large u behaviors of the integrands of the integrals that occur in the 
density density correlation function one can examine the behaviors of solutions to the mode 
equation in this limit. It turns out that these can be obtained by using a Volterra series to 



find solutions for the mode functions in the L and R regions. To solve Eq. (6.2) in terms of 
a Volterra series one first notes that Ves vanishes in the limit x — )■ 0. If Xo is any solution 
to the equation when Vcs = then a formal solution to the full equation is 



1 



X 



X^(x) = x^(x) - - / dy* sm[uix* - y*)]Vcsiy*)x''iy*) ■ (7-5) 

^ J-oo 

This equation can be solved by iteration. After the first iteration one finds 

X'^ix) = xoi^) - - r dy* sm[uix* - y*)]V,siy*)Xoiy*) ■ (7-6) 



In the out region Ves vanishes in the limit x ^ oo. In this case if Xo is a solution when 
Ves = then the formal solution is 

X^(x) = Xoi^) + - dy* sm[u{x* - y*)]Kfr(l/*)xf , (7-7) 



After one iteration 



X^(x) = xoi^) + - I dy* sm[co{x* - y*)]V,f,{y*)x^{y*) . (T.i 



00 



X* 



To obtain the large u limit of the mode functions we note that if Ve// = then 

Xc = cos(u;x*) , (7.9a) 

= sin(cux*) , (7.9b) 

= cos{ux*) , (7.9c) 

xf = sin(wx*) . (7.9d) 



Substituting these expressions into Eqs. (7.6) and (7.8), and using simple trigonometric 
identities one finds that to 0(1/0;) 

(7.10a) 
(7.10b) 
(7.10c) 
(7.10d) 

Jx 
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= COS 


[ux*) - 


sm{ujx*) 


r dyV{y), 
Jo 


2uj 




= sini 


ux*) + 


cosiux*) 


r dyV{y), 
Jo 


2uj 


xf 


= cos 


[OJX*) + 


sin(cjx*) 


POD 

/ dyViy), 

J X 


2uj 


Xs 


= sini 


UJX*) - 


COs{(jJX*) 


POO 

/ dyViy). 

J X 


2u 



One can use these results to obtain the large u limits of the matching parameters in 



Eqs. (6.5). However from Eqs. (6.8) it can be seen that what one really needs are the sums 
and differences of A and -D, and of B and C . To 0{uj^) one finds that A = D = 1 and 
B = C = Q. To Oiuj^"^) the sums and differences are 



A + D 
A-D 



2 




B+C = , 
B-C = - 

CO 



dyV{y) . 



(7.11) 



These results can be substituted into the integrals for the density density correlation 
function in terms of the modes x ^^id the matching parameters A, B, C, and D to obtain 
an approximate expression which is valid in the large u limit. This is then subtracted from 
the density density correlation function and added back on. The part that is added back on 
is computed analytically by first computing the approximate two point function and then 
taking the relevant derivatives. In the rest one has integrals which can now be computed 
numerically because they are well behaved in the limit a; — > oo. 

VIII. NUMERICAL RESULTS 



Numerical computations of the density-density correlation function have been carried out 
for the sound speed profile 



c= \icl + -{cl-cl) 



' 2 ^ . fx + b 
1 + -tan-i 



with 



b = ay tan 



The values used for some of the constants were 



3 

4 ' 
1 

2 ' 

C2 = 1 . 



^0 
Cl 



U) 



(8.2) 



(8.3) 
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This is the same type of profile as that used in Ref. [10] where numerical calculations of 
the density density correlation function were carried out in the context of condensed matter 
physics. However for those calculations Ci = 1 and C2 = 1/2 so that the BEC was moving 
to the right, not the left [20]. Also 6 = so that the horizon was not at x = 0. 

Calculations were done for several values of the parameter c^^, ranging from 1/4 to 8. 
Typical values used for the infrared cutoff A in the integrals over the frequency u ranged 
from 2 X 10~^ to 2 x 10^^. All of the plots shown here have an infrared cutoff of u; = 2 x 10~^. 



A careful analysis of (8.1) shows that the derivatives of c(x) are larger near the horizon 



for smaller values of cr^. The hydrodynamic approximation that we are using is valid only 



if the derivatives of c(x) are relatively small. For the profile in (8.1) the hydrodynamical 
approximation was shown in [10] to be valid when ^ 4. Thus while we have computed 
the density density correlation function for smaller values of (jy we display our results only 
for cr^ = 8. 

In all cases the correlation function was computed for equal values of the lab time T and in 



fact specifically for T = 0. The arbitrary constants in Eq. (3.3) and in Eq. (3.10) for the out 
region have been given the values Xi = X3 = 1. In the in region the corresponding constants 
were given the values X2 = = —1. To have continuity of the coordinate v = t + x* across 
the future horizon it is easy to show by equating the expressions for v in the in and out 
regions at a; = that 

--Tt^. (8.4) 



As can be seen in Eqs. (5.6) and (5.9), the two point function can be written in terms 
of an integral (/) over the mode functions and and an integral (J) over the mode 
functions uf^^. Thus the density density correlation function can also be separated into 
terms containing integrals over and and terms containing integrals over -u^"*. 

In Ref. [9] the density density correlation function was computed analytically with the 
following assumptions and restrictions: It was computed with one point inside the horizon 
and one point outside. The assumption was made that the point inside was in a region 
where the sound speed was constant and the point outside was in a region where the sound 
speed was also constant but had a different value. Having neglected the effective potential 



Veff in Eq. (3.12) and hence the backscattering, the computed correlation corresponds to 
that in Fig. [9] (between the modes in Fig. |6]and Fig. [?]). The main result of that paper is 
the existence of a negative correlation peak (i.e. a trough) in the correlation function which 
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would not be there if there was no horizon. 

In Ref. [TU] a fully quantum mechanical computation of the density density correlation 
function was done and a comparison was made with the result of [S]. The existence of 
the negative correlation peak when one point is inside and one point is outside the horizon 
was confirmed. The peak value of the correlation function along with the full width at half 
maximum of the peak were compared when the points were far enough from the horizon that 
the sound speed was effectively constant. They were found to be in approximate agreement 
only for > 4. Also found were two other correlation peaks, both of which are substantially 
weaker than the one found in P|. One of these correlation peaks is negative (between the 
modes in Fig. [7]and Fig. [s] modes, see Fig. 11). It occurs when both points are inside the 
horizon. The other is a very weak (positive) correlation peak which occurs when one point 
is outside and one point is inside the horizon (between the modes in Fig. |6] and Fig. Isl see 



Fig. 10). 



Our results for the density density correlation function when ct„ = 8 are shown in Fig. 12 



When both points are inside or both points are outside the horizon there is a divergence 
in the correlation function when the points come together. Also the correlation function is 
symmetric about the point separation. Thus only half the range of values of the points is 
shown in those plots. There is no divergence in the density density correlation function and 
it is not symmetric about the point separation when one point is inside and the other point 
is outside the horizon. Thus we show the full range of values of the points in this case. 

When both points are outside the horizon the only feature we observed in the density 
density correlation function is the divergence as the two points come together. 

When one point is inside and one point is outside the horizon we find the same negative 
correlation peak as was found when the potential is zero. It it interesting in this case to 
look separately at the contribution of the modes and those of the uf modes. These 



are shown in the first two plots in Fig. 13 Note that the contribution from the uf modes is 
actually a positive correlation peak. However, it is significantly smaller in magnitude than 
the negative correlation peak from the other set of modes. The effect on the magnitude of 
the negative correlation peak by including the effective potential in the mode equations and 
including the contribution from the uf modes was observed to be small, usually less than 
10%. 



As can be seen in the bottom plot of Fig. 13 there is also a small positive correlation peak 
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FIG. 12: The density density correlation function is shown for = 8. Both points are outside 
the horizon in the top left plot, one point is inside and one point is outside the horizon in the top 
right plot, and both points are inside the horizon in the bottom plot. In all cases the lab time T is 
the same for both points and and an infrared cutoff in the frequency integrals of A = 2 x 10"'^ has 
been used. Note that in the top left and bottom plots the bottom looks flat. This is a distortion 
caused by imposing a limit on the range of values of the correlation function that are plotted. 



which occurs when one point is inside and the other point is outside the horizon. This one 
lies closer to the horizon than the larger negative one and its maximum value is significantly 
smaller than that of the main peak. This appears to correspond to the correlator between 
the modes in Fig. [6] and Fig. [8} see Fig. 10 A careful examination of Fig. 13 shows that the 
contribution from the mJ§ and modes is a positive correlation peak and the contribution 
from the m^"* modes is a negative correlation peak. However, the two only partially cancel 
when added together and the result is the small positive correlation peak in the bottom plot. 



In the bottom plot of Fig. 12 one finds a negative correlation peak which occurs when 
both points are inside the horizon. This corresponds to the correlator between the modes in 
Fig. [7] and Fig. |8} see Fig. 11 If the vertical scale is decreased by a factor of ten in order to 



show more detail then as shown in Fig. 14 the density density correlation function reaches 
a maximum value on either side of the negative correlation peak. The outer maximum 
appears to be due to the existence of the negative correlation peak just mentioned and also 
the very large one that occurs when the points come together. There must be a maximum 
in between. Something similar may happen for the other maximum which is in between the 
negative correlation peak and the horizon. However, we are not able to numerically compute 
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FIG. 13: The density density correlation function is shown for = 8 when one point is inside 
and one point is outside the horizon. The plot on the top left shows the contribution from the 
u^^ and n^* modes. The plot on the top right shows the contribution from the uj^^ modes. The 
plot on the bottom shows the total. In all cases the lab time T is the same for both points and 
and an infrared cutoff in the frequency integrals of A = 2 x 10~^ has been used. The vertical scale 
is set so that the positive correlation peak on the upper left plot, the two maxima on the upper 
right plot, and the positive correlation peak in the bottom plot can be seen. This has the effect 
of distorting the shape of the negative correlation peak in the plots on the top left and bottom 
so that it appears flat near its maximum. Note that the vertical scale of the bottom plot is 1/10 
as large as the top plots. It is necessary to do this to clearly show the small peak because of the 
partial cancelation in this peak that occurs when the two contributions are added together. 

the density density correlation function when one point is on the horizon so we don't know 
its value there. 

When one point is inside and one point is outside the horizon, a quantitative comparison 
of our numerical results with the analytic calculation in [9] gives agreement to within approx- 
imately 10% or better for the size of the negative correlation peak for all values of that 
were considered. As mentioned previously the comparison of the numerical computations 
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done in [TO] with the analytic calculation in [9] showed approximate agreement for cr„ ^ 3. 
For smaller values the size of the peak predicted by the analytic calculation becomes much 
larger than that found in the numerical calculations in ^U\. Since our numerical results 
take the effective potential Ves in the mode equation into account, they give the correct 
quantum field theory in curved space prediction for the density density correlation function. 
Our numerical calculations thus confirm that the discrepancy is due to the fact that the 
hydrodynamical approximation breaks down for < 3 due to the rapid variation of the 
sound speed near the horizon for those values. 

From the above discussion it is clear that our numerical results reproduce all of the same 
late time peaks in the correlation function that were found in [ini UHl [Ej. This is true even 
for (Tt, ^ 3. In [in] only the results for the negative correlation peak when one point is inside 
and one point is outside the horizon were shown for more than one value of a^. For the other 
peaks the results shown in [TOl [131 E] were only shown for the case = 1/2. However, we 
have been given access [20] to some of the numerical data that was generated for [TU] . This 
has allowed quantitative comparisons to be made between our numerical results and those 
numerical results for the other correlation peaks. 

In the numerical data we were given [20j we could only unambiguously identify the positive 
correlation peak when one point is inside and one point is outside the horizon for cr^ < 4 
although there was evidence for it for larger values of a^. In those cases our value for the 
size of the peak was always somewhat larger. However, as discussed below the effects of the 
infrared cutoff are large enough that the size of this peak can be significantly affected. 

For the negative correlation peak which occurs when both points are inside the horizon 
we find that at cr^ > 4 there is agreement to roughly 25% or better between the two data 
sets with that agreement being better for larger values of as expected. For all the values 
of ay considered our results give a larger value for the size of the negative correlation peak 
than those for |10j . 



As mentioned above, and as seen in Fig. [14| when both points are inside the event horizon 
we find evidence for a maximum in the density density correlation function on either side 
of the negative correlation peak. In the numerical data we were given [20j the maximum 
that is farther from the horizon than the negative peak is always seen. There is also always 
evidence for the maximum we see that is closer to the horizon but for cr^, > 2 an actual 
maximum is not seen. Instead the value of the correlation function just keeps increasing as 
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FIG. 14: The density density correlation function is shown for = 8 when both points are inside 



the horizon. This is the same as the bottom plot in Fig. 12 but the vertical scale has been changed 
to more clearly show the maxima on either side of the negative correlation peak. This has the 
effect of distorting the shape of the larger negative correlation peak so that it appears flat near its 



maximum. 



one point approaches the horizon (and the other point remains far from the horizon). For 
(T„ = 1 a maximum is seen. In some cases for our data the inner maximum has a positive 
value and is thus a positive correlation peak. The outer maximum however has a negative 
value except for = 1 where it has a very small positive value. In the other data set the 
outer maximum is less than zero for cr^ < 6 and greater than zero for larger values of a^. For 
(Tt, > 1 the outer maximum is smaller than the inner maximum. Neither maximum was seen 
in [in] but both of them can be seen in [13] , although only the larger, inner one is identified 
as a correlation peak. 

As discussed previously the infrared cutoff that we use does have an effect on the value of 
the density density correlation function and will have a strong effect for values of the cutoff 
that are either too large or too small. Different values of the infrared cutoff were considered 
for some of our calculations. It was found for example for = 8 that when a comparison 
was made between calculations with a cutoff of A = 2 x 10~^ and a calculation with a cutoff 
of A = 2 X 10^^ that the differences were typically in the range (2 — 4) x 10^®. This is 
small enough that the negative correlation peaks are not significantly affected. However, 
the positive correlation peak and the two maxima discussed above are close enough to zero 
that their values can be affected. A similar comparison between the cutoffs of A = 2 x 10^^ 
and A = 2 X 10""^ yields differences that typically range from (1 — 4) x 10~^ which is about 
an order of magnitude higher. Thus a cutoff of A = 2 x 10"'^ is probably too high. 
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IX. CONCLUSIONS 



We have developed a method to numerically compute the density density correlation 
function for a BEC which serves as an analogue black hole. In the process we have shown 
that it is necessary to incorporate an infrared cutoff in the frequency. So long as this cutoff 
is small in the sense that its product with the point separation is small, then it does not 
significantly affect the main peaks or troughs in the correlation function. However it does 
significantly affect the value of the correlation function in regions where the magnitude of 
that function is small. 

In order to correctly compute the correlation function it is necessary to consider modes of 
arbitrarily high frequencies. This is exactly what one expects when working with an effective 
field theory [21] . Unless this is done ultraviolet cutoff effects dominate the structure of the 
correlation function. 

To take these modes into account it was found that a type of regularization must be 
used in order to compute the correlation function numerically. One must subtract off terms 
whose amplitude grows or stays the same as the frequency increases. Then these terms are 
added back again and the integral over u is computed formally by starting with an integral 
of the form 



and taking derivatives with respect to u. This reason for this is related to the fact that for the 
original correlation function one must first compute the integral over the mode frequencies 
and then take the derivatives necessary to compute the density density correlation function. 
These two operations do not commute. But if one performs a regularization as just described 
then there is no problem in first computing the derivatives and then integrating over the 
frequencies. This is the desired order of operations for numerical computations. 

Our numerical results show that for the models considered the size and location of the 
main negative peak in the correlation function when one point is inside and one point is 
outside the horizon is not changed significantly by including the potential. However, includ- 
ing the potential does correctly reproduce the features which were found in the condensed 
matter calculations in jlOj, namely a negative correlation peak which occurs when both 




(9.1) 



or 




(9.2) 
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points are inside the horizon and a small positive correlation peak that occurs when one 
point is inside and the other point is outside the horizon. Finally we find two maxima which 
occur on either side of the negative correlation peak when both points are inside the horizon. 
These maxima do not appear on the plots of Ref. [lOj, but they do appear in the plot in [13] 
although only the one closer to the horizon is identified as a correlation peak. 

In cases where the hydrodynamical approximation is valid there is reasonably good quan- 
titative agreement with the calculations described in Ref. [10] which used sophisticated nu- 
merical condensed matter simulations. The late time features found in those calculations 
can be reproduced using just the gravity analogy and quantum field theory in curved space 
techniques. 
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